Growth - Midway Chondria tumulosa

Data collection was divided into Canopy and Understory

Dictionary of terms:

plant_part = canopy or understory (under),

nitrate = nutrient dose, which includes phosphate,

d9_growth_percent = calculated growth from final weight-initial weight divided by initial,

mean_mins28 = mean of summed time in minutes for 9 day period each bin in canopy or understory experienced temperatures 28C or higher

Midway Experiment

Prepare the data

Load packages

library(lme4)
library(lmerTest)
library(effects)
library(car)
library(MuMIn)
library (dplyr)
library(emmeans)
library(DHARMa)
library(performance)
library(patchwork)
library(rstatix)
#for plots and tables
library(ggplot2)
library(ggpubr)
library(forcats)
library(RColorBrewer)
library(tidyverse)
library(sjPlot)
library(sjmisc)
library(gt)
library(purrr)
library(stringr)
library(tidyr)
library(piecewiseSEM)
library(easystats)
library(magrittr)

Open weight dataset

mid_growth <- read.csv("../../../data/midway_2024/input/midway_growth.csv")

Calculate growth rate from initial and final weights, make new columns

mid_growth$d9_growth_percent <- round(((mid_growth$final_weight - mid_growth$initial_weight) / mid_growth$initial_weight * 100), digits = 2)
mid_growth$d5_growth_percent <- round(((mid_growth$d5_weight - mid_growth$initial_weight) / mid_growth$initial_weight * 100), digits = 2)

Add a new column that gets rid of characters in ID column

mid_growth$plant_ID <- as.factor(substr(mid_growth$ID, 4, 5))

Make a new column for run using the third value in the ID, then change the 1 to 2 and 2 to 3

mid_growth$run <- substr(mid_growth$ID, 3, 3)
mid_growth$run <- as.factor(ifelse(mid_growth$run == 1, 2, 3)) #actual run numbers are 2 and 3

Assign treatment as characters from integers then to factors

mid_growth$nitrate <- as.factor(as.character(mid_growth$treatment))

Get rid of the ppt in salinity

mid_growth$salinity <- substr(mid_growth$salinity, 1, 2)

Combine nitrate and salinity for a treatment number

mid_growth <- mid_growth %>%
  mutate(treatment = case_when(
    nitrate == 1 & salinity == 35 ~ "1",
    nitrate == 1 & salinity == 28 ~ "2",
    nitrate == 2 & salinity == 35 ~ "3",
    nitrate == 2 & salinity == 28 ~ "4",
    nitrate == 3 & salinity == 35 ~ "5",
    nitrate == 3 & salinity == 28 ~ "6",
    nitrate == 4 & salinity == 35 ~ "7",
    nitrate == 4 & salinity == 28 ~ "8",
  ))

Create subsets for the plots

Use Day 9 for final analysis

canopy_g <- subset(mid_growth, plant_part == "canopy" & d9_growth_percent != 42.45)
canopy_g$treatment_graph[canopy_g$treatment == 1] <- "1) 0.5umol/35 ppt"
canopy_g$treatment_graph[canopy_g$treatment == 2] <- "2) 0.5umol/28 ppt"
canopy_g$treatment_graph[canopy_g$treatment == 3] <- "3) 2umol/35 ppt" 
canopy_g$treatment_graph[canopy_g$treatment == 4] <- "4) 2umol/28 ppt"
canopy_g$treatment_graph[canopy_g$treatment == 5] <- "5) 4umol/35 ppt"
canopy_g$treatment_graph[canopy_g$treatment == 6] <- "6) 4umol/28 ppt"
canopy_g$treatment_graph[canopy_g$treatment == 7] <- "7) 8umol/35 ppt"
canopy_g$treatment_graph[canopy_g$treatment == 8] <- "8) 8umol/28 ppt"

under_g <- subset(mid_growth, plant_part == "under" & d9_growth_percent != -62.82 & 
                    d9_growth_percent != -32.50) #removing one influential datapoint
under_g$treatment_graph[under_g$treatment == 1] <- "1) 0.5umol/35 ppt"
under_g$treatment_graph[under_g$treatment == 2] <- "2) 0.5umol/28 ppt"
under_g$treatment_graph[under_g$treatment == 3] <- "3) 2umol/35 ppt" 
under_g$treatment_graph[under_g$treatment == 4] <- "4) 2umol/28 ppt"
under_g$treatment_graph[under_g$treatment == 5] <- "5) 4umol/35 ppt"
under_g$treatment_graph[under_g$treatment == 6] <- "6) 4umol/28 ppt"
under_g$treatment_graph[under_g$treatment == 7] <- "7) 8umol/35 ppt"
under_g$treatment_graph[under_g$treatment == 8] <- "8) 8umol/28 ppt"

Add new column to subsets from “time over 28” summary dataset

time_over_28 <- read_csv("../../../data/midway_2024/transformed/mins_28_plant_part.csv") #load dataset

Keep only relevant columns, convert to factors before combining datasets so data types match

mean_mins28 <- time_over_28 %>%
  select(plant_part, run, mean_mins28_plant_part) %>% #keep only relevant columns of data
  mutate(run = as.factor(run), mean_mins28_plant_part = as.factor(mean_mins28_plant_part))

Join the time over 28 dataset to the growth dataset using run and plant_part as matches CANOPY

canopy_g <- canopy_g %>%
  left_join(mean_mins28, by = c("run", "plant_part")) #join the datasets

canopy_g <- canopy_g %>%
  rename(mean_mins28 = mean_mins28_plant_part) #name is too long
glimpse(canopy_g)

## Rows: 96
## Columns: 14
## $ ID                <chr> "ct101", "ct102", "ct103", "ct104", "ct105", "ct106"…
## $ treatment         <chr> "1", "2", "3", "4", "5", "6", "7", "8", "1", "2", "3…
## $ salinity          <chr> "35", "28", "35", "28", "35", "28", "35", "28", "35"…
## $ plant_part        <chr> "canopy", "canopy", "canopy", "canopy", "canopy", "c…
## $ initial_weight    <dbl> 0.3964, 0.3709, 0.4576, 0.4222, 0.3401, 0.3468, 0.42…
## $ d5_weight         <dbl> 0.4140, 0.3949, 0.4280, 0.4321, 0.3892, 0.3963, 0.43…
## $ final_weight      <dbl> 0.4159, 0.4040, 0.4143, 0.4307, 0.4150, 0.4028, 0.47…
## $ d9_growth_percent <dbl> 4.92, 8.92, -9.46, 2.01, 22.02, 16.15, 12.14, -42.45…
## $ d5_growth_percent <dbl> 4.44, 6.47, -6.47, 2.34, 14.44, 14.27, 2.29, -10.94,…
## $ plant_ID          <fct> 01, 02, 03, 04, 05, 06, 07, 08, 09, 10, 11, 12, 13, …
## $ run               <fct> 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2…
## $ nitrate           <fct> 1, 1, 2, 2, 3, 3, 4, 4, 1, 1, 2, 2, 3, 3, 4, 4, 1, 1…
## $ treatment_graph   <chr> "1) 0.5umol/35 ppt", "2) 0.5umol/28 ppt", "3) 2umol/…
## $ mean_mins28       <fct> 2825, 2825, 2825, 2825, 2825, 2825, 2825, 2825, 2825…

Join the time over 28 dataset to the growth dataset using run and plant_part as matches UNDERSTORY

under_g <- under_g %>%
  left_join(mean_mins28, by = c("run", "plant_part")) #join the datasets

under_g <- under_g %>%
  rename(mean_mins28 = mean_mins28_plant_part) #name is too long
glimpse(under_g)

## Rows: 94
## Columns: 14
## $ ID                <chr> "ct149", "ct150", "ct151", "ct152", "ct153", "ct154"…
## $ treatment         <chr> "1", "2", "3", "4", "5", "6", "7", "8", "1", "2", "3…
## $ salinity          <chr> "35", "28", "35", "28", "35", "28", "35", "28", "35"…
## $ plant_part        <chr> "under", "under", "under", "under", "under", "under"…
## $ initial_weight    <dbl> 0.4615, 0.4468, 0.3672, 0.4865, 0.4462, 0.4084, 0.48…
## $ d5_weight         <dbl> 0.4983, 0.4398, 0.3940, 0.4496, 0.4498, 0.4422, 0.46…
## $ final_weight      <dbl> 0.5332, 0.4092, 0.3755, 0.4294, 0.4265, 0.4649, 0.40…
## $ d9_growth_percent <dbl> 15.54, -8.42, 2.26, -11.74, -4.42, 13.83, -16.43, 0.…
## $ d5_growth_percent <dbl> 7.97, -1.57, 7.30, -7.58, 0.81, 8.28, -3.42, 2.08, 0…
## $ plant_ID          <fct> 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 61, 62, …
## $ run               <fct> 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2…
## $ nitrate           <fct> 1, 1, 2, 2, 3, 3, 4, 4, 1, 1, 2, 3, 3, 4, 4, 1, 1, 2…
## $ treatment_graph   <chr> "1) 0.5umol/35 ppt", "2) 0.5umol/28 ppt", "3) 2umol/…
## $ mean_mins28       <fct> 1398, 1398, 1398, 1398, 1398, 1398, 1398, 1398, 1398…

Analysis – Canopy

Run linear mixed-effects model and run checks for fit

Make histograms

canopy_g %>% 
ggplot(aes(d9_growth_percent)) + #day 9
  geom_histogram(binwidth=5, fill = "goldenrod1", color = "black", linewidth = 0.25, alpha = 0.85) +
  theme_bw()

canopy_g %>% ggplot(aes(d5_growth_percent)) + #day 5
  geom_histogram(binwidth=5, fill = "goldenrod1", color = "black", linewidth = 0.25, alpha = 0.85) +
  theme_bw()

Run the model

Mean_min28 used as random effect in lieu of run, but neither had an effect, also rlc_order had no effect

growth_model_canopy <- lmer(formula = d9_growth_percent ~ salinity + nitrate +
                            (1 | plant_ID), data = canopy_g, REML = TRUE)

Check model fit

isSingular(growth_model_canopy)

## [1] FALSE

hist(resid(growth_model_canopy))

plot(resid(growth_model_canopy) ~ fitted(growth_model_canopy))

qqnorm(resid(growth_model_canopy))
qqline(resid(growth_model_canopy))

Check the performance of the model for dataset: canopy

performance::check_model(growth_model_canopy)

r.squaredGLMM(growth_model_canopy)

##             R2m       R2c
## [1,] 0.09262105 0.3014886

summary(growth_model_canopy)

## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: d9_growth_percent ~ salinity + nitrate + (1 | plant_ID)
##    Data: canopy_g
## 
## REML criterion at convergence: 738.7
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -3.03933 -0.63855  0.03797  0.56284  2.36387 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  plant_ID (Intercept)  39.51    6.286  
##  Residual             132.12   11.495  
## Number of obs: 96, groups:  plant_ID, 48
## 
## Fixed effects:
##             Estimate Std. Error     df t value Pr(>|t|)  
## (Intercept)    5.639      3.316 43.000   1.701   0.0963 .
## salinity35     5.960      2.966 43.000   2.009   0.0508 .
## nitrate2      -3.383      4.195 43.000  -0.807   0.4243  
## nitrate3       2.340      4.195 43.000   0.558   0.5798  
## nitrate4      -5.145      4.195 43.000  -1.227   0.2266  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##            (Intr) slnt35 nitrt2 nitrt3
## salinity35 -0.447                     
## nitrate2   -0.632  0.000              
## nitrate3   -0.632  0.000  0.500       
## nitrate4   -0.632  0.000  0.500  0.500

View random effects levels

tab_model(growth_model_canopy, show.intercept = TRUE, show.se = TRUE, show.stat = TRUE, show.df = TRUE, show.zeroinf = TRUE)

	d9_growth_percent
Predictors	Estimates	std. Error	CI	Statistic	p	df
(Intercept)	5.64	3.32	-0.95 – 12.23	1.70	0.093	89.00
salinity [35]	5.96	2.97	0.07 – 11.85	2.01	0.048	89.00
nitrate [2]	-3.38	4.19	-11.72 – 4.95	-0.81	0.422	89.00
nitrate [3]	2.34	4.19	-5.99 – 10.68	0.56	0.578	89.00
nitrate [4]	-5.15	4.19	-13.48 – 3.19	-1.23	0.223	89.00
Random Effects
σ2	132.12
τ00 plant_ID	39.51
ICC	0.23
N plant_ID	48
Observations	96
Marginal R2 / Conditional R2	0.093 / 0.301

plot(allEffects(growth_model_canopy))

Construct null model to perform likelihood ratio test (REML must be FALSE)

canopy_nitrate_null <- lmer(formula = d9_growth_percent ~ salinity + 
                              (1 | plant_ID), data = canopy_g, REML = FALSE)
canopy_model2 <- lmer(formula = d9_growth_percent ~ nitrate + salinity +
                        (1 | plant_ID), data = canopy_g, REML = FALSE)
anova(canopy_nitrate_null, canopy_model2)

## Data: canopy_g
## Models:
## canopy_nitrate_null: d9_growth_percent ~ salinity + (1 | plant_ID)
## canopy_model2: d9_growth_percent ~ nitrate + salinity + (1 | plant_ID)
##                     npar    AIC    BIC  logLik deviance  Chisq Df Pr(>Chisq)
## canopy_nitrate_null    4 770.61 780.86 -381.30   762.61                     
## canopy_model2          7 772.50 790.45 -379.25   758.50 4.1109  3     0.2497

canopy_salinity_null <- lmer(formula = d9_growth_percent ~ nitrate + 
                               (1 | plant_ID), data = canopy_g, REML = FALSE)
canopy_model3 <- lmer(formula = d9_growth_percent ~ nitrate + salinity + 
                        (1 | plant_ID), data = canopy_g, REML = FALSE)
anova(canopy_salinity_null, canopy_model3)

## Data: canopy_g
## Models:
## canopy_salinity_null: d9_growth_percent ~ nitrate + (1 | plant_ID)
## canopy_model3: d9_growth_percent ~ nitrate + salinity + (1 | plant_ID)
##                      npar   AIC    BIC  logLik deviance Chisq Df Pr(>Chisq)  
## canopy_salinity_null    6 774.8 790.19 -381.40    762.8                      
## canopy_model3           7 772.5 790.45 -379.25    758.5 4.308  1    0.03793 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Plot canopy growth

canopy_growth_plot <- canopy_g %>% 
  ggplot(aes(treatment_graph, d9_growth_percent)) + 
  geom_boxplot(size=0.5) + 
  geom_point(alpha = 0.75, size = 5, aes(color = salinity), position = "jitter", show.legend = TRUE) + 
  labs(x="Nitrate", y= "9-Day Growth (%)", title= "A", subtitle = "Chondria tumulosa -- Canopy") + 
  scale_x_discrete(labels = c("0.5 μmol N", "0.5 μmol N", "2 μmol N", "2 μmol N", 
                              "4 μmol N", "4 μmol N", "8 μmol N", "8 μmol N")) + 
  ylim(-45, 50) + stat_mean() + 
  scale_color_manual(values = c("goldenrod1", "darkgoldenrod3")) +
  geom_hline(yintercept=0, color = "red", linewidth = 0.5, alpha = 0.5) +
  theme_bw() +
  theme(legend.position.inside = c(0.90,0.90), plot.title = element_text(face = "bold", vjust = -15, hjust = 0.05), 
        plot.subtitle = element_text(face = "italic", size = 14, vjust = -20, hjust = 0.05))
canopy_growth_plot

Analysis – Understory

Make a histogram of the growth rate data

under_g %>% ggplot(aes(d9_growth_percent)) + #day 9
  geom_histogram(binwidth=5, fill = "#990066", color = "black", linewidth = 0.25, alpha = 0.85) +
  theme_bw()

under_g %>% ggplot(aes(d5_growth_percent)) + #day 5
  geom_histogram(binwidth=5, fill = "#990066", color = "black", linewidth = 0.25, alpha = 0.85) +
  theme_bw()

Mean_min28 used as random effect in lieu of run, RLC order had no effect

growth_model_under <- lmer(formula = d9_growth_percent ~ salinity + nitrate +
                              (1 | plant_ID) + (1 | mean_mins28), data = under_g, REML = TRUE)

isSingular(growth_model_under)

## [1] FALSE

hist(resid(growth_model_under))

plot(resid(growth_model_under) ~ fitted(growth_model_under))

qqnorm(resid(growth_model_under))
qqline(resid(growth_model_under))

Check the performance of the model for dataset: under

performance::check_model(growth_model_under)

r.squaredGLMM(growth_model_under)

##             R2m       R2c
## [1,] 0.07896071 0.2748286

summary(growth_model_under)

## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: d9_growth_percent ~ salinity + nitrate + (1 | plant_ID) + (1 |  
##     mean_mins28)
##    Data: under_g
## 
## REML criterion at convergence: 653.7
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -2.50427 -0.51736 -0.02634  0.54235  2.21577 
## 
## Random effects:
##  Groups      Name        Variance Std.Dev.
##  plant_ID    (Intercept)  0.9344  0.9666  
##  mean_mins28 (Intercept) 18.7592  4.3312  
##  Residual                72.9123  8.5389  
## Number of obs: 94, groups:  plant_ID, 48; mean_mins28, 2
## 
## Fixed effects:
##             Estimate Std. Error     df t value Pr(>|t|)   
## (Intercept)   13.207      3.646  1.704   3.623  0.08643 . 
## salinity35    -1.844      1.785 42.742  -1.033  0.30730   
## nitrate2      -7.434      2.524 42.720  -2.946  0.00519 **
## nitrate3      -3.589      2.524 42.720  -1.422  0.16218   
## nitrate4      -3.287      2.496 41.754  -1.317  0.19510   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##            (Intr) slnt35 nitrt2 nitrt3
## salinity35 -0.245                     
## nitrate2   -0.335 -0.015              
## nitrate3   -0.335 -0.015  0.490       
## nitrate4   -0.342  0.000  0.495  0.495

View random effects levels

tab_model(growth_model_under, show.intercept = TRUE, show.se = TRUE, show.stat = TRUE, show.df = TRUE, show.zeroinf = TRUE)

	d9_growth_percent
Predictors	Estimates	std. Error	CI	Statistic	p	df
(Intercept)	13.21	3.65	5.96 – 20.46	3.62	<0.001	86.00
salinity [35]	-1.84	1.78	-5.39 – 1.70	-1.03	0.304	86.00
nitrate [2]	-7.43	2.52	-12.45 – -2.42	-2.95	0.004	86.00
nitrate [3]	-3.59	2.52	-8.61 – 1.43	-1.42	0.159	86.00
nitrate [4]	-3.29	2.50	-8.25 – 1.68	-1.32	0.191	86.00
Random Effects
σ2	72.91
τ00 plant_ID	0.93
τ00 mean_mins28	18.76
ICC	0.21
N plant_ID	48
N mean_mins28	2
Observations	94
Marginal R2 / Conditional R2	0.079 / 0.275

plot(allEffects(growth_model_under))

Construct null model to perform likelihood ratio test REML must be FALSE

under_nitrate_null <- lmer(formula = d9_growth_percent ~ salinity + 
                              (1 | plant_ID) + (1 | mean_mins28), data = under_g, REML = FALSE)
under_model2 <- lmer(formula = d9_growth_percent ~ nitrate + salinity +
                        (1 | plant_ID) + (1 | mean_mins28), data = under_g, REML = FALSE)
anova(under_nitrate_null, under_model2)

## Data: under_g
## Models:
## under_nitrate_null: d9_growth_percent ~ salinity + (1 | plant_ID) + (1 | mean_mins28)
## under_model2: d9_growth_percent ~ nitrate + salinity + (1 | plant_ID) + (1 | mean_mins28)
##                    npar    AIC    BIC  logLik deviance Chisq Df Pr(>Chisq)  
## under_nitrate_null    5 689.38 702.10 -339.69   679.38                      
## under_model2          8 686.74 707.08 -335.37   670.74 8.644  3    0.03442 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

under_salinity_null <- lmer(formula = d9_growth_percent ~ nitrate + (1 | plant_ID) + (1 | mean_mins28), data = under_g, REML = FALSE)
under_model3 <- lmer(formula = d9_growth_percent ~ nitrate + salinity + (1 | plant_ID) + (1 | mean_mins28), data = under_g, REML = FALSE)
anova(under_salinity_null, under_model3)

## Data: under_g
## Models:
## under_salinity_null: d9_growth_percent ~ nitrate + (1 | plant_ID) + (1 | mean_mins28)
## under_model3: d9_growth_percent ~ nitrate + salinity + (1 | plant_ID) + (1 | mean_mins28)
##                     npar    AIC    BIC  logLik deviance  Chisq Df Pr(>Chisq)
## under_salinity_null    7 685.88 703.68 -335.94   671.88                     
## under_model3           8 686.74 707.08 -335.37   670.74 1.1413  1     0.2854

Plot understory growth

under_growth_plot <- under_g %>% 
  ggplot(aes(treatment_graph, d9_growth_percent)) + 
  geom_boxplot(size=0.5) + 
  geom_point(alpha = 0.75, size = 5, aes(color = salinity), position = "jitter", show.legend = TRUE) + 
  labs(x="Nitrate", y= "9-Day Growth (%)", title= "B", subtitle = "Chondria tumulosa -- Understory") + 
  scale_x_discrete(labels = c("0.5 μmol N", "0.5 μmol N", "2 μmol N", "2 μmol N", 
                              "4 μmol N", "4 μmol N", "8 μmol N", "8 μmol N")) + 
  ylim(-35, 40) + stat_mean() + 
  scale_color_manual(values = c("maroon", "maroon2")) +
  geom_hline(yintercept=0, color = "red", linewidth = 0.5, alpha = 0.5) +
  theme_bw() +
  theme(legend.position.inside = c(0.90,0.90), plot.title = element_text(face = "bold", vjust = -15, hjust = 0.05), 
        plot.subtitle = element_text(face = "italic", size = 14, vjust = -20, hjust = 0.05))
under_growth_plot

Summarize the means for nitrate and salinity

mid_growth %>% group_by(nitrate, plant_part) %>% summarise_at(vars(d9_growth_percent), list(mean = mean))

## # A tibble: 8 × 3
## # Groups:   nitrate [4]
##   nitrate plant_part  mean
##   <fct>   <chr>      <dbl>
## 1 1       canopy      8.62
## 2 1       under      12.3 
## 3 2       canopy      5.24
## 4 2       under       2.12
## 5 3       canopy     11.0 
## 6 3       under       7.07
## 7 4       canopy      3.47
## 8 4       under       9.00

mid_growth %>% group_by(salinity, plant_part) %>% summarise_at(vars(d9_growth_percent), list(mean = mean))

## # A tibble: 4 × 3
## # Groups:   salinity [2]
##   salinity plant_part  mean
##   <chr>    <chr>      <dbl>
## 1 28       canopy      4.09
## 2 28       under       7.45
## 3 35       canopy     10.1 
## 4 35       under       7.79